This is my script to try and analyse data from scratch - DATASET 2 -> MERSCOPE
library(Seurat)
Loading required package: SeuratObject
Loading required package: sp
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
(status 2 uses the sf package in place of rgdal)
Attaching package: ‘SeuratObject’
The following objects are masked from ‘package:base’:
intersect, t
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
library(ggplot2)
library(scCustomize)
scCustomize v2.1.2
If you find the scCustomize useful please cite.
See 'samuel-marsh.github.io/scCustomize/articles/FAQ.html' for citation info.
library(readr)
library(pheatmap)
library(matrixStats)
library(spdep)
Loading required package: spData
To access larger datasets in this package, install the spDataLarge package with:
`install.packages('spDataLarge', repos='https://nowosad.github.io/drat/', type='source')`
Loading required package: sf
Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
library(geojsonR)
ImageDimPlot(seurat,fov='COLON')
ImageDimPlot(seurat,fov='COLON')
# Mark the blank probes that are detected in each cell - don't think we have the other codeword categories that are present
seurat[["Negative.Control.Codeword"]] <- CreateAssayObject(counts =data$transcripts[grepl('Blank',rownames(data$transcripts)),])
ImageFeaturePlot(seurat, "nCount_MERSCOPE", axes = T) + scale_fill_viridis_c()
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
# crop only does image, not dataset
#cropped <- Crop(seurat[["COLON"]], x = c(8000, 11000), y = c(2500, 5500), coords = "plot")
#seurat[["ROIIW"]] <- cropped
# MERSCOPE CELL GETS read in as an integer but the integers are too big - need to read in cell ids as character so use tidyverse to read in set as character
global_coordinates <- data.frame(seurat@images$COLON$centroids)
global_coordinates
## This does not work - need to order correctly
seurat$global_X <- global_coordinates$x
seurat$global_Y <- global_coordinates$y
seurat_subset <- seurat[,seurat$global_X < 11000 & seurat$global_X > 8000 & seurat$global_Y < 5500 & seurat$global_Y > 2500]
Warning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objects
ImageFeaturePlot(seurat_subset, "nCount_MERSCOPE", axes = T) + scale_fill_viridis_c()
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
ImageFeaturePlot(seurat_subset, "nFeature_MERSCOPE") + scale_fill_viridis_c()
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
ImageFeaturePlot(seurat_subset, "volume") + scale_fill_viridis_c()
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
seurat_subset[["SIZE_FILTER_LARGE"]] <- seurat_subset$volume < quantile(seurat_subset$volume, .99)
ImageDimPlot(seurat_subset, group.by="SIZE_FILTER_LARGE")
seurat_subset[["SIZE_FILTER_SMALL"]] <- seurat_subset$volume > quantile(seurat_subset$volume, .01)
ImageDimPlot(seurat_subset, group.by="SIZE_FILTER_SMALL")
```{r}
Error: attempt to use zero-length variable name
seurat_subset$TRANSCRIPT_FILTER <- seurat_subset$nCount_MERSCOPE >= 15
ImageDimPlot(seurat_subset, group.by="TRANSCRIPT_FILTER")
ImageFeaturePlot(seurat_subset, "nCount_Negative.Control.Codeword") + scale_fill_viridis_c()
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
ImageFeaturePlot(seurat_subset, "nCount_Negative.Control.Probe") + scale_fill_viridis_c()
Error in `FetchData()`:
! None of the requested variables were found: nCount_Negative.Control.Probe
Backtrace:
1. Seurat::ImageFeaturePlot(seurat_subset, "nCount_Negative.Control.Probe")
3. SeuratObject:::FetchData.Seurat(object = object, vars = c(features, split.by[1L]), cells = cells)
ImageFeaturePlot(seurat_subset, "nFeature_Negative.Control.Codeword") + scale_fill_viridis_c()
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
seurat_subset_filtered <- subset(seurat_subset, SIZE_FILTER_LARGE & SIZE_FILTER_SMALL & TRANSCRIPT_FILTER)
Warning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objects
seurat_subset_filtered <- SCTransform(seurat_subset_filtered, assay = "MERSCOPE", clip.range = c(-10, 10))
Running SCTransform on assay: MERSCOPE
Running SCTransform on layer: counts
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 550 by 91161
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 549 genes, 5000 cells
Found 63 outliers - those will be ignored in fitting/regularization step
Second step: Get residuals using fitted parameters for 550 genes
Computing corrected count matrix for 550 genes
Calculating gene attributes
Wall clock passed: Time difference of 15.35141 secs
Determine variable features
Centering data matrix
|
| | 0%
|
|===============================================================================================| 100%
Getting residuals for block 1(of 19) for counts dataset
Getting residuals for block 2(of 19) for counts dataset
Getting residuals for block 3(of 19) for counts dataset
Getting residuals for block 4(of 19) for counts dataset
Getting residuals for block 5(of 19) for counts dataset
Getting residuals for block 6(of 19) for counts dataset
Getting residuals for block 7(of 19) for counts dataset
Getting residuals for block 8(of 19) for counts dataset
Getting residuals for block 9(of 19) for counts dataset
Getting residuals for block 10(of 19) for counts dataset
Getting residuals for block 11(of 19) for counts dataset
Getting residuals for block 12(of 19) for counts dataset
Getting residuals for block 13(of 19) for counts dataset
Getting residuals for block 14(of 19) for counts dataset
Getting residuals for block 15(of 19) for counts dataset
Getting residuals for block 16(of 19) for counts dataset
Getting residuals for block 17(of 19) for counts dataset
Getting residuals for block 18(of 19) for counts dataset
Getting residuals for block 19(of 19) for counts dataset
Centering data matrix
|
| | 0%
|
|===============================================================================================| 100%
Finished calculating residuals for counts
Set default assay to SCT
seurat_subset_filtered <- RunPCA(seurat_subset_filtered)
PC_ 1
Positive: EPCAM, LGR5, EPHB3, RGMB, SOX9, CDH1, PKM, MYC, VEGFA, MUC1
AXIN2, CDCA7, IDH1, MKI67, ERBB2, PROX1, ASCL2, ERBB3, SMOC2, CTNNB1
EPHB4, NOTCH1, PCNA, HDAC1, EPHA2, LGALS9, MET, JUN, MCM2, CEACAM1
Negative: COL1A1, ACTA2, FN1, COL4A1, COL5A1, PDGFRB, PECAM1, MYH11, CD248, ETS1
CAV1, PLVAP, ITGA1, HLA-DRA, ITGA5, DES, ENG, MMP11, ZEB1, ANGPT2
DUSP1, WWTR1, AKT3, ELN, COL11A1, NRP1, ADAMTS4, SFRP2, VWF, CSF1
PC_ 2
Positive: SPP1, HLA-DRA, ITGB2, CYBB, FCGR3A, C1QC, CD14, CSF1R, HLA-DRB1, MRC1
HLA-DQA1, LYZ, PTPRC, CIITA, HLA-DPB1, ITGAX, CXCR4, FCGR2A, HLA-DMA, SERPINA1
ITGAM, HAVCR2, MAFB, IL4I1, TREM2, CD4, TLR2, DUSP1, CEBPB, FOS
Negative: FN1, COL4A1, ACTA2, COL1A1, SMOC2, PDGFRB, COL5A1, LGR5, RGMB, ITGA1
CD248, MMP11, PROX1, CAV1, EPHB3, ANGPT2, CDCA7, MYH11, PLVAP, NOTCH1
ELN, ADAMTS4, CTNNB1, INSR, ETS1, NDUFA4L2, MKI67, ASCL2, WWTR1, SOX9
PC_ 3
Positive: SMOC2, MKI67, LGR5, CDCA7, PLK1, RGMB, PCNA, PROX1, BIRC5, HLA-DRA
CCNB1, FOXM1, MCM6, MCM2, FN1, CDK4, EPHB3, ITGB2, C1QC, CTNNB1
CSF1R, HLA-DRB1, MYBL2, CYBB, ASCL2, HLA-DQA1, FCGR3A, AURKB, EZH2, CD14
Negative: FOS, VEGFA, EPHA2, CEACAM1, JUNB, EGR1, LAMC2, SLC26A3, CDKN1A, MUC1
KIT, LAMB3, CXCL1, LRP1, ATF3, NFKBIA, EPCAM, LDHA, PLOD2, PPARD
PKM, IFNGR2, COL4A1, DUSP1, NFKB2, INSR, TCF7L2, KDR, MMP7, ITGB1
PC_ 4
Positive: FN1, LGR5, SMOC2, RGMB, ERBB3, PROX1, COL1A1, MMP11, LRP1, COL5A1
STAT6, ERBB2, COL11A1, DES, RORC, EPHB3, SFRP2, CDKN1B, SOX9, AXIN2
EPCAM, TNFSF10, LRP5, ELN, VEGFA, JUN, BCL2, IGF1R, FGFR3, EPHA4
Negative: MKI67, PLK1, BIRC5, PCNA, FOXM1, PLVAP, CCNB1, PECAM1, PKM, KDR
MYBL2, VWF, MMRN2, AURKB, ENG, MCM6, AURKA, MCM2, FLT4, COL4A1
FLT1, ANGPT2, EZH2, PDGFB, ADAMTS4, NRP1, ETS1, BRCA1, CLEC14A, E2F1
PC_ 5
Positive: PLVAP, PECAM1, INSR, VWF, MMRN2, KDR, ENG, LGR5, FLT4, SMOC2
FLT1, CDH5, PDGFB, ETS1, RGMB, ANGPT2, NRP1, CLEC14A, ADAMTS4, PREX2
COL4A1, ERBB3, PROX1, CD40, CXCR4, ACKR3, NOTCH1, ITGA5, THBD, TNFRSF4
Negative: MKI67, COL1A1, COL5A1, PLK1, ACTA2, FN1, PKM, BIRC5, PCNA, FOXM1
DES, CCNB1, EGR1, FOS, COL11A1, TGFBI, MYH11, SFRP2, TNC, ELN
CDKN1A, MYBL2, AURKB, AURKA, MCM6, MCM2, JUNB, CSF1, CA9, MMP11
ElbowPlot(seurat_subset_filtered, 50)
PC_Plotting(seurat_subset_filtered, dim_number = 8)
FeaturePlot(seurat_subset_filtered, "MRC1", reduction = "pca") + scale_color_viridis_c()
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
ImageFeaturePlot(seurat_subset_filtered, "PC_1") + scale_fill_viridis_c()
Warning: No FOV associated with assay 'SCT', using global default FOVScale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
ImageFeaturePlot(seurat_subset_filtered, "MRC1", size=.5) + scale_fill_viridis_c()
Warning: No FOV associated with assay 'SCT', using global default FOVScale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
#Add annotation
ref <- readRDS("/project/shared/spatial_data_camp/datasets/SINGLE_CELL_REFERENCES/COLON_HC_5K_CELLS.RDS")
ref
DimPlot(ref)
ref <- SCTransform(ref, residual.features =rownames(seurat_subset_filtered_UMAP))
ref <- RunPCA(ref)
ref <- RunUMAP(ref, dims=1:20)
DimPlot(ref, label=T, repel=T)
ps <- AggregateExpression(ref, features = rownames(seurat_subset_filtered_UMAP), normalization.method = "LogNormalize", assays="RNA", return.seurat = T)
#Transfer labels
ImageDimPlot(seurat_subset_filtered_UMAP, size=.5)
Warning: No FOV associated with assay 'SCT', using global default FOV
FeaturePlot(seurat_subset_filtered_UMAP, "CD3E", label=T, repel=T)+ scale_color_viridis_c(direction=-1)
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
top
[1] "LGR5" "RORC" "RGMB" "SMOC2" "ERBB3" "MKI67" "PLK1"
[8] "AURKB" "CCNB1" "FOXM1" "SLC26A3" "MMP7" "CXCL1" "CEACAM1"
[15] "LAMC2" "FN1" "CDH5" "PROX1" "BCL2" "PDGFRB" "COL5A1"
[22] "COL1A1" "ELN" "DES" "CLCA1" "MUC2" "SERPINA1" "Blank-46"
[29] "PTGDR2" "SPP1" "CSF1R" "HLA-DRA" "ITGB2" "MMP12" "PLVAP"
[36] "VWF" "PECAM1" "MMRN2" "FLT4" "CD3E" "TRAC" "CD2"
[43] "CCL5" "EOMES"
Clustered_DotPlot(seurat_subset_filtered_UMAP, features = top, group.by = "predicted.id", k=9)
[[1]]
[[2]]
ImageDimPlot(seurat_subset_filtered_UMAP, fov="ROI", boundaries="segmentation", border.color = "black" )
Error in ImageDimPlot(seurat_subset_filtered_UMAP, fov = "ROI", boundaries = "segmentation", :
No compatible spatial coordinates present
ImageDimPlot(seurat_subset_filtered_UMAP, group.by = "predicted.id")
Warning: No FOV associated with assay 'SCT', using global default FOV
cell_colours <- c("#F8766D", "#DB8E00", "#AEA200", "#64B200", "#00BD5C", "#00C1A7",
"#00BADE", "#00A6FF", "#B385FF", "#EF67EB", "#FF63B6")
names(cell_colours) <- c("Epithelium", "Fibroblasts", "T-Cells", "Myofibroblasts", "Macrophages",
"Glia", "Endothelium", "Telocytes", "Plasma", "B-Cells", "Pericytes")
cell_colours
Epithelium Fibroblasts T-Cells Myofibroblasts Macrophages
"#F8766D" "#DB8E00" "#AEA200" "#64B200" "#00BD5C"
Glia Endothelium Telocytes Plasma B-Cells
"#00C1A7" "#00BADE" "#00A6FF" "#B385FF" "#EF67EB"
Pericytes
"#FF63B6"
ImageFeaturePlot(seurat_subset_filtered_UMAP, "SMOC2") + scale_fill_viridis_c()
Warning: No FOV associated with assay 'SCT', using global default FOVScale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
markers <- FindAllMarkers(seurat_subset_filtered_UMAP, group.by="predicted.id")